The evolution of environmentally mediated social interactions and posthumous spite under isolation by distance

Many social interactions happen indirectly via modifications of the environment, e.g. through the secretion of functional compounds or the depletion of renewable resources. Here, we derive the selection gradient on a quantitative trait affecting dynamical environmental variables that feed back on reproduction and survival in a finite patch-structured population subject to isolation by distance. Our analysis shows that the selection gradient depends on how a focal individual influences the fitness of all future individuals in the population through modifications of the environmental variables they experience, weighted by the neutral relatedness between recipients and the focal. The evolutionarily relevant trait-driven environmental modifications are formalized as the extended phenotypic effects of an individual, quantifying how a trait change in an individual in the present affects the environmental variables in all patches at all future times. When the trait affects reproduction and survival through a payoff function, the selection gradient can be expressed in terms of extended phenotypic effects weighted by scaled relatedness. We show how to compute extended phenotypic effects, relatedness, and scaled relatedness using Fourier analysis, which allow us to investigate a broad class of environmentally mediated social interactions in a tractable way. We use our approach to study the evolution of a trait controlling the costly production of some lasting commons (e.g. a common-pool resource or a toxic compound) that can diffuse in space and persist in time. We show that indiscriminate posthumous spite readily evolves in this scenario. More generally, whether selection favours environmentally mediated altruism or spite is determined by the spatial correlation between an individual’s lineage and the commons originating from its patch. The sign of this correlation depends on interactions between dispersal patterns and the commons’ renewal dynamics. More broadly, we suggest that selection can favour a wide range of social behaviours when these have carry-over effects in space and time.


Introduction
Organisms continually interact with one another in ways that significantly impact their survival and reproduction.Such social interactions are incredibly diverse.Still, they can usefully be classified as to whether they occur directly between individuals, such as grooming or fighting over resources, or as to whether they are indirectly mediated by environmental modifications, such as through the depletion or enrichment of resources, or the release or detoxification of pollutants [1].Direct social interactions thus take place among contemporaries who are physically close to one another, while environmentally mediated interactions can extend much further in space and time.When environmental modifications have long-lasting and long-ranging effects, indirect social interactions may occur between individuals whose lifetimes show little or no overlap.This may lead to forms of trans-generational helping (e.g. when the underconsumption of a slowly renewable resource ensures healthy stock maintenance for future generations) or harming (e.g. when overconsumption leads to stock collapse and poor harvest in the future).
An extreme form of trans-generational social behaviour is posthumous altruism, which is a behaviour that results in a net reduction of the survival and/or reproduction of an actor to benefit only recipients living beyond the actor's death.Several models have explored the biological scenarios that favour the evolution of posthumous altruism [2,3].But perhaps more striking is posthumous spite, which is a behaviour that results in a net reduction of the survival and/or reproduction of an actor to harm recipients that live only after the actor's death.The conditions that may favour posthumous spite remain unclear.Understanding this, and more generally the evolution of environmentally mediated social behaviour, requires describing how selection shapes traits that underlie both direct and indirect social interactions.
The theory devoted to the evolution of quantitative traits influencing direct social interactions is well established (see, e.g.[4][5][6][7][8][9]).One of the main contributions of this body of work has been to highlight the importance of limited dispersal for determining how selection shapes social traits in populations that are subdivided into finite groups and spatially structured [4,5].Under limited dispersal, stochastic demographic effects resulting from finite group (or patch or deme) size generate genetic associations, whereby individuals expressing the same traits may be more or less likely to interact directly with one another than with individuals expressing alternative traits.The importance of such associations is enshrined in the fact that the selection gradient on a quantitative trait can be expressed as a marginal (or gradient) form of Hamilton's rule [4,5].This gradient, which captures the first-order effects of selection, is sufficient to determine the trait values towards which a population converges under mutation-limited evolution ( [7], i.e. to characterise convergence stability [10]).
In contrast, modelling the evolution of social interactions mediated by abiotic or biotic environmental variables is significantly more challenging in spatially structured populations.This is because computing the selection gradient in this case also requires computing the joint distribution of pairwise relatedness and environmental variables in the population (under neutrality [5,25]).Generally, this distribution is the stationary solution to a high-dimensional stochastic dynamical system that is difficult to analyse or, in some cases, even to characterise.The challenge is apparent in models that allow for trait-driven changes in local demography.Even in the island model of dispersal, where the spatial structure is only implicit [26], there is typically no analytical solution to the distribution of demographic states within groups.As a result, evaluating the selection gradient on traits with effects on local demography and/or ecology is computationally challenging (e.g.[25,[27][28][29][30]).This "curse of dimensionality" becomes even more acute under isolation by distance, as the size of the state space on which relevant environmental (or demographic) variables fluctuate blows up exponentially, with the selection gradient now requiring the distribution of states among as well as within groups (e.g.eq.22 in [25]).
To circumvent this challenge, two approximations have been suggested.One is the pair approximation developed for lattice-structured populations, where typically at most one individual lives in sites connected by stepping-stone dispersal [31][32][33][34][35][36][37][38].Pair approximation is based on moment equations of the demographic state distribution that ignore third and higher-order moments.Under this approximation, the selection gradient can be written in the form of Hamilton's marginal rule, thus allowing for a sharp understanding of some of the effects of demography on the evolution of social behaviour ( [19,20] see also [39]).However, pair approximation is not straightforward when considering arbitrarily complex dispersal patterns (e.g.[40]), patches with more than one individual, or trait-driven environmental state variables.
Another approximation relies on considering the dynamics of environmental state variables to be deterministic with a stable fixed point, so that there are no environmental stochastic fluctuations in the absence of genetic variation [41].The selection gradient can then be readily expressed as a marginal Hamilton's rule with inter-temporal fitness effects arising through trait-driven environmental modifications at different temporal distances.In addition to being simpler to compute than the original problem, this decomposition allows the delineation between a component of selection resulting from direct social interactions and a component arising indirectly through changes in the environmental dynamics.So far, this approach has been applied only to the island model-hence, in the absence of isolation by distance [41,42].For populations exhibiting isolation by distance, there exist general, marginal Hamilton's rule like, formulas.These express the selection gradient in terms of inter-temporal and spatial effects of trait expression on the fitness of all possible recipients [2,43].However, how environmental modifications mediate the long-lasting and long-ranging fitness effects due to trait expression remains implicit in these general formulas.To better understand selection resulting from indirect social interactions via environmental feedback, fitness effects need to be unpacked in terms of trait-driven environmental modifications at different temporal and spatial distances.
Here, we do exactly that.We fully characterise the selection gradient on a trait that impacts the deterministic dynamics of environmental state variables that can be abiotic or biotic, and which feed back on the survival and reproduction of the evolving species in a general model of isolation by distance.Using Fourier analysis, we express this gradient in terms of extended phenotypic effects and relatedness coefficients scaled to local competition, both of which provide biological insights about the nature of selection and are straightforward to compute for a wide range of classical evolutionary models (e.g. the Wright-Fisher model and the Cannings model).We use our results to investigate the evolution of environmentally mediated helping and harming through space and time, such as via the production of a lasting common-pool resource or the release of a durable toxic compound.Our analyses indicate that indiscriminate spite where individuals suffer a cost to harm others living in the future, even beyond the actor's death, readily evolves by natural selection.This result is in contrast with previous findings suggesting that indiscriminate spite rarely evolves when it occurs through direct interactions.More broadly, our model and results suggest that selection can favour a wide range of social behaviours when these are mediated in space and time through environmental feedbacks.

Spatial structure, life cycle, traits and environmental variables
We consider a population of homogeneous individuals subdivided among D homogeneous patches (or demes or groups), each carrying N adult individuals (see Table 1 for a list of the key symbols used).The population is censused at discrete demographic time steps between which the following events occur in cyclic order: (a) reproduction and adult survival; (b) dispersal among patches; and (c) density-dependent regulation within patches such that each patch contains N adult individuals at the beginning of the next demographic time step.
Patches are arranged homogeneously in d dimensions, with D j patches in dimension j 2 {1, . .., d}.For example, under a lattice structure in a one-dimensional habitat, D = D 1 patches are arranged on a circle, while in a two-dimensional habitat, D = D 1 × D 2 patches are arranged on a torus.We denote by G ¼ fði 1 ; i 2 ; . . .; i d Þ : 0 � i j < D j g the set of all patches.The fact that patches are arranged homogeneously means that, at a technical level, we can endow the set G with an abelian group structure (see e.g.[44] for the use of such group structure in evolutionary biology), which will allow us to leverage techniques from Fourier analysis on finite abelian groups (Box 1).
Each patch is characterized by a quantitative state variable representing a biotic or abiotic environmental factor, which we refer to as an environmental state variable (e.g., density of a common-pool resource or pollutant, habitat quality).Meanwhile, each individual in the population is characterised by a genetically determined quantitative trait (e.g.consumption of a resource, release of a pollutant, investment into habitat maintenance) that influences the environment and the individual's survival and reproduction.We are interested in the evolution of this trait under the following three assumptions.
1. Trait and environmentally mediated reproduction and survival.By expressing the evolving trait, individuals can directly affect the survival and/or reproduction of any other individual in the population.For example, individuals may engage in costly fights for resources in other patches and return to their own patch to share these resources with patch neighbours.
The effects of trait expression on others are assumed to be: (a) spatially invariant, i.e. the marginal effect of an individual from patch i ¼ ði 1 ; i 2 ; . . .; i d Þ 2 G on the survival and/or reproduction of an individual in patch j ¼ ðj 1 ; j 2 ; . . .; j d Þ 2 G only depends on the "distance" j − i between the two patches (where j − i is calculated from the abelian group operation, see Box 1); and (b) spatially symmetric, i.e. the marginal effect of an individual from patch i on an individual in patch j is equal to the effect from j to i.We refer to these two characteristics (a) and (b) together as spatial homogeneity.The survival and reproduction of an individual may also depend on the environmental state variable of each patch,

Box 1. Fourier analysis on finite abelian groups
We assume that the set of patches G is endowed with an abelian group structure, which allows us to consider more general spatial structures than just lattice models (e.g.hierarchical structures).The group G is defined as the direct product, Here, we followed the convention of population genetics (e.g.[5,46,47]) and defined the Fourier transform in terms of the character χ k (h) (instead of its conjugate, given in Eq I.E, which is more standard in mathematics and engineering).As such, the Fourier transform of f gives the characteristic function of f when f is a probability distribution.For instance, the Fourier transform MðhÞ ¼ P k2G m k w k ðhÞ is the characteristic function of the dispersal distribution m k .The original function is recovered by using F ðhÞ� w k ðhÞ; ðI:DÞ where L k ðF Þ is the inverse transform of F at k 2 G, which is defined in terms of the conjugate of χ k (h): ! ðI:EÞ (e.g.[82]).Another property that we use in our analysis is the orthogonality relation between characters, i.e., the identity (see [82], p. 169).also in a spatially homogeneous way (i.e. the marginal effect of the environmental state variable of a patch i on the survival and reproduction of an individual residing in patch j only depends on the distance j − i, and is equal to the effect from j to i).

2.
Dispersal.Each individual either stays in its natal patch or disperses to another patch.Dispersal occurs with non-zero probability so that patches are never completely isolated.We assume that dispersal is spatially homogeneous, i.e. the probability of dispersal from one patch i to another j depends only on the distance j − i = k between the two patches (spatial invariance), and is equal to the probability of dispersing the distance i − j = −k (spatial symmetry).We can thus write m k = m −k for the probability that an individual disperses to a patch at a distance k from its natal patch (with P k2G m k ¼ 1).

Trait and environmentally mediated environmental dynamics. Through trait expression,
individuals can affect environmental state variables from one demographic time step to the next.For example, the environmental variable may be a common-pool resource that individuals absorb locally, or a pollutant produced by individuals which then diffuses in the environment.Such trait effects on the environment are also spatially homogeneous (i.e. the marginal effect of an individual from patch i on the environmental state variable of patch j only depends on the distance j − i and is equal to the effect from j to i).These trait-driven environmental modifications can thus lead to inter-temporal, environmentally mediated social interactions.

The focal individual, its fitness, and environmental dynamics
The spatial homogeneity underlying all processes described above means that the patch indexed as 0 2 G can be taken as a representative patch, and that any individual in this patch can be taken as a representative individual from the population.We refer to this patch and to this individual as, respectively, the focal patch and the focal individual.In general, we refer to patch k 2 G at time (or "generation") t � 0 prior to the focal generation as patch "k, t".Thus, the focal patch corresponds to patch 0, 0. In the following, we introduce some notation to describe trait and environmental variation in the population (see Fig 1C is the expected number of successful offspring (i.e.offspring that establish as adults, including the surviving self) produced over one demographic time by the focal individual with trait z • , when the trait average among other individuals at the different spatial positions is z 0,0 , and environmental state variables are n 0,0 .These state variables are obtained from the solution to the system of equations .The map g depends on (i) the traits in the whole population expressed at the previous generation via z k,t+1 (recall that t goes back in time), which is a circular permutation of the elements of z 0,t+1 with z k,t+1 as first element [e.g. for a one dimensional lattice, z 0,t+1 = (z 0,t+1 , z 1,t+1 , . .., z D−1,t+1 ), z 1,t+1 = (z 1,t+1 , z 2,t+1 , . .., z 0,t+1 ), z 2,t+1 = (z 2,t+1 , z 3,t+1 , . .., z 0,t+1 , z 1,t+1 ), and so on]; and (ii) the environmental state variables of all patches at the previous generation via n k,t+1 .Due to the recursive nature of Eq (2), the environmental state variables in the focal generation, n 0,0 , depend on the whole history of traits z H = (z 0,1 , z 0,2 , . ..) expressed in the population prior to the focal generation.As a result, the fitness of a focal individual may also depend on the traits expressed by all other previous individuals across space and time.To make this dependence explicit, we write the fitness of the focal individual as w(z • , z 0,0 , n 0,0 (z H )). We make the additional assumption that in a monomorphic population where all individuals express the same resident phenotype z, the deterministic environmental dynamics described by the map g have a unique hyperbolically stable equilibrium point, identical in each patch, and satisfying where z = (z, . .., z) and n ¼ ðn; . . .; nÞ are vectors of dimension D whose entries are all equal to trait value z and environmental state variable n, respectively.This is sometimes called the spatially homogeneous or flat solution in multi-patch ecological systems (p.235 in [45]).The existence of such a solution entails that, in the absence of genetic variation, all patches converge to the same environmental equilibrium n, which may depend on the resident trait z.
One useful property of a monomorphic population at such an equilibrium is that fitness must be equal to one, i.e. wðz; z; nÞ ¼ 1 holds.This is because the total population size is constant, and consequently, each individual exactly replaces itself on average.The fitness function ( 1) and the recursion (2) assume that fitness and the environmental dynamics can be written as functions of trait averages within patches.This said, this assumption does not limit us to only considering situations where effects within patches are additive.Indeed, because we are interested in convergence stability and thus in the first-order effects of selection (i.e. the first-order effects of trait variation), Eqs (1) and ( 2) are sufficient to model biological scenarios with arbitrary complicated non-additive traits effects among individuals within and/or between patches, for instance through complementarity or antagonism.Just a little care may be required when defining these expressions from an individual-based model (p.95 in [5]).

Evolutionary dynamics
We assume that the quantitative trait evolves through rare mutations of small phenotypic effects, such that the evolutionary dynamics proceeds as a trait substitution sequence on the state space Z � R (i.e. the process of "long-term evolution" for finite populations described in [7]).We are interested in characterising convergence stable trait values, which are local attractors of the trait substitution sequence.To do so, we base ourselves on the first-order effects of selection on the fixation probability of a mutant that arises as a single copy in a population monomorphic for a resident trait value [7,11,12].Technical details about our derivations can be found in S1 Text and accompanying boxes.Our main findings are summarized below.

Recipient-centered perspective: Intra-and inter-temporal fitness effects
holds, where the function referred to as the selection gradient, can be written as the sum of two terms, given by and (see Appendix A in S1 Text for a derivation).Here, all quantities are evaluated in a monomorphic population with all individuals expressing the resident trait value z, and the environmental state variable in all patches are at the environmental equilibrium n (Eq 3).A trait value z* satisfying condition (4) constitutes a candidate end-point of evolution.More specifically, z* is a mode of the stationary phenotypic distribution under the trait substitution sequence (see eq. A-5 in S1 Text and e.g.[7] for details).It is also possible for a mode of this distribution to lie on the boundary of a closed trait space, e.g. of Z ¼ ½z min ; z max � so that a mode sits at z min and/or z max .For a mode to be at z min , s(z min ) < 0 must hold, in which case z min is convergence stable.Similarly, if s(z max ) > 0 holds, then z max is convergence stable and a mode of the stationary phenotypic distribution.Both s w (z) and s e (z) depend on (i) marginal fitness effects, i.e. on derivatives of focal fitness (that we interpret below), and (ii) the relatedness R k,t between the focal individual and another randomly sampled individual from patch k, t.Such relatedness coefficient is defined as where μ is the mutation rate at the evolving locus; Q k,t is the stationary probability that an allele sampled in the focal individual is identical by descent with a homologous allele sampled in another individual chosen at random from patch k, t under neutrality (i.e. in a population monomorphic for z); and � Q t ¼ P k2G Q k;t =D is the average probability of identity by descent between two homologous alleles sampled in two individuals living t generations apart.The probability of identity by descent Q k,t , and thus the relatedness coefficient R k,t , may depend on the resident phenotype z, but we leave this dependence implicit for readability.
Relatedness R k,t quantifies the extent to which an individual that is sampled from patch k, t is more (when R k,t > 0) or less (when R k,t < 0) likely than a randomly sampled individual to carry an allele identical by descent to one carried by the focal individual at a homologous locus.To illustrate this notion, consider a Wright-Fisher process (where there is no adult survival and individuals are semelparous), which is the reference model for probabilities of identity by descent under isolation by distance (e.g.[5,46,47]).For this model, the relatedness coefficients R k,t for t = 1, 2, 3, . . .can be shown to be given by where and MðhÞ is the Fourier transform (or characteristic function) of the m k probabilities (see Box 1 for definitions of Fourier transforms, character functions χ k (h), and their inverses � w k ðhÞ; Appendix B in S1 Text for an example of the characteristic function of a dispersal distribution; and [48] for a derivation of Eq 8).The relatedness coefficient between two individuals in the same generation, R k,0 , is given by Eq (8) with t = 2 (i.e.R k,0 = R k,2 holds for all k).In a panmictic or randomly mixing population (where m k = 1/D for all k), relatedness between any two individuals is zero (i.e.R k,t = 0 for all k and all t; as the Fourier transform reduces to MðhÞ ¼ 1 if h = 0 and 0 otherwise using property I.F in Box 1).However, as soon as dispersal is non-uniform (i.e. if m k 6 ¼ m j for some k 6 ¼ j), relatedness varies among individuals from different patches according to spatial and temporal distance.In particular, when dispersal is limited so that individuals have a tendency to remain in their natal patch, relatedness between the focal individual and individuals in the same patch from the same generation increases (R 0,0 > 0).Because the average relatedness is zero (i.e.P k2G R k;t =D ¼ 0 holds from Eq 7), the focal individual must also be negatively related to individuals residing in at least one other patch (i.e.R k,0 < 0 must hold for some k 6 ¼ 0).Which patches those are depends on patterns of dispersal.Under short-range dispersal, the focal individual tends to be positively related to individuals in patches nearby and negatively related to individuals further away (Figs 2C and 3C).Under long-range dispersal, relatedness can be negative between individuals living in patches at intermediate distances (Figs 2D and 3D).In Eq (6), the derivative @w/@z • is the effect of a trait change in the focal individual on its own fitness.Likewise, @w/@z k,t is the effect of a trait change in the whole set of individuals living in patch k, t on the fitness of the focal individual.Such effect is weighted by relatedness R k,t in the expressions given by Eq (6).As such, s w (z) (Eq 6a) is the net effect of all intra-temporal effects on fitness, while s e (z) (Eq 6b) is the net effect of all inter-temporal effects (i.e.all effects within and between demographic periods, respectively).More broadly, Eq (6) consists in the sum of effects on the fitness of a focal individual stemming from all individuals that currently live (Eq 6a) or have lived (Eq 6b) in the population.Here, the focal individual is the recipient of phenotypic effects, present and past (recipient-centered perspective).How past phenotypic effects are mediated by environmental dynamics is left implicit in Eq (6), contained in Eq (6b) through n 0,0 (z H ). In the next section, we expose such environmental effects by unpacking Eq (6b) and shifting from a recipient-centered to an actor-centered perspective.

Actor-centered perspective: Environmentally mediated extended phenotypic effects
To understand natural selection on social traits, it is often helpful to see the focal individual as the actor, rather than the recipient, of phenotypic effects [4,5,49].To shift to this perspective, we can leverage the space-time homogeneity of our model to see that @w/@z k,t in Eq ( 6) is equivalent to the total effect of the focal individual on the fitness of the individuals in a patch at distance k at t steps in the future, and that relatedness R k,t (Eq 8) quantifies the extent to which an individual sampled in a patch at distance k at t steps in the future is more (or less) likely than a randomly sampled individual to carry an allele identical by descent to one in the focal individual at a homologous locus [2].These considerations readily lead to an actor-centered perspective for selection on intra-temporal effects, s w (z) (Eq 6a).
For selection on inter-temporal effects, s e (z), we further need to unpack the phenotypic effects through the environmental dynamics in Eq (6b).To do so, we now let the time index t � 0 denote time forward so that n k,t is the value of the environmental state variable in patch k at t time steps in the future of the focal generation (t = 0), and likewise let z k,t denote the collection of population phenotypes at time t in the future.Environmental dynamics forward in time are characterised by rewriting Eq (2) as where z R k;0 is equal to z k,0 except that the component z 0,0 in this vector is replaced with z R 0;0 ¼ z � =N þ ðN À 1Þz 0;0 =N, i.e. the average phenotype in the focal patch including the focal individual (e.g. for a one-dimensional lattice, z R k;0 ¼ ðz k;0 ; z kþ1;0 ; . . .; z DÀ 1;0 ; z R 0;0 ; z 1;0 ; . . .; z kÀ 1;0 Þ).Eq (9) brings upfront the potential complexity of characterising the environmental consequences of a trait change in the focal individual.This is because the trait of the focal individual, z • , influences the environmental state variable of potentially any patch k over one generation, n k,1 , which can in turn have knock-on effects in the future on n k,2 , n k,3 , and so on throughout space in an interactive way.To characterise such effects, we write for the extended phenotypic effect of the focal individual on the environmental state variable in patch k at t generations in the future.Selection on inter-temporal effects s e (z) (Eq 6b) can then be written in terms of these extended phenotypic effects as s e ðzÞ ¼ X (see Appendix C.1 in S1 Text for a derivation).Eq (11) indicates that s e (z) consists in the total effect of the focal individual on the fitness of individuals in each patch k, t in the future, via a change in the environmental state of possibly all patches j − k, t, where fitness is weighted by their relatedness R k,t to the focal individual.From the point of view of the focal individual, relatedness R k,t can then be thought of as the "genetic value" of an individual randomly sampled in patch k, t in units of fitness.More specifically, R k,t can be interpreted as the number of units of its own fitness that the focal individual is willing to exchange with an individual from patch k, t against one unit of theirs without changing the mutant's probability of fixation at z*.
Selection thus favours the focal individual sacrificing some of its own fitness to increase fitness in patch k, t when R k,t > 0, and to decrease fitness when R k,t < 0. How such sacrifice impacts the environment encountered by recipients is quantified by the extended phenotypic effect e j−k,t in Eq (11).
The remaining challenge is how to compute e k,t , given the complex repercussions that a perturbation has in time and space (i.e.how to quantify a perturbation in the coupled dynamical system defined by Eq 9).We show in Appendix C.2 in S1 Text that this can be achieved through Fourier analysis using the following building blocks.First, we let ψ k be the focal individual's effect on the environmental state variable of patch k over one generation.Owing to our space-time homogeneity assumptions, this effect can be calculated as In practice, the rightmost expression is often more useful than the expression between equal signs in Eq (12) (as the rightmost expression only requires characterising the environmental map gðz R 0;0 ; n 0;0 Þ of the focal patch, e.g.Eq 26 below).Similarly, we let be the effect of the environmental state variable of one patch on the environmental state variable of another patch at distance k over one generation.With the above notation, and writing C(h) and CðhÞ for the Fourier transforms of ψ k and c k respectively, the extended phenotypic effect can be efficiently computed as the inverse Fourier transform of CðhÞ tÀ 1 CðhÞ ¼ E t ðhÞ (see Box 1 for the definition of inverse Fourier transforms), namely |ffl ffl ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl ffl ffl } The form of E t ðhÞ indicates that e k,t can be considered a perturbation in the dynamics of an environmental state variable that ripples into the future.This perturbation has its origin in a focal individual whose trait affects the environmental state variables of possibly multiple patches over one generation (captured by C(h) in Eq 14).This one-generational change then propagates through space over t − 1 generations owing to the environmental dynamics, finally impacting the environment of individuals t generations downstream of the focal individual (captured by CðhÞ tÀ 1 in Eq 14).In Box 2, we generalize Eqs ( 11)-( 14) to multi-dimensional environmental dynamics, that is, when multiple environmental state variables can be affected by the evolving trait and whose dynamics can interact with one another (e.g., metacommunity dynamics).Eq (11) together with Eqs ( 12)-( 14) constitutes a basis to quantify and understand selection on traits that have inter-temporal effects through the environment under isolation by distance.These equations formalise the intuition behind inclusive fitness arguments for environmentally mediated social interactions, saying that natural selection tends to favor traits whose environmental effects benefit the fitness of relatives (i.e.individuals more likely to carry genes that are identical by descent).Here, the fundamental currency is individual fitness, and its exchange rate among individuals is given by relatedness.However, in many cases it is not fitness that is directly impacted by traits or the environment, but rather some intermediate payoff, such as energy, the amount of prey caught, or the size of a breeding territory.In turn, this payoff influences survival and reproduction, which together determine fitness.We next explore such scenarios.

Payoff-mediated fitness: Scaled relatedness or the genetic value of others in payoff units
3.3.1 Payoff and fitness.Following much of evolutionary game theory (see e.g.[8] for a textbook treatment), we now consider cases where fitness depends on a payoff function summarizing social interactions (e.g. a case where the energy that an individual accrues depends on its foraging behaviour, the foraging behaviour of others, and how resources are distributed in the environment).We let this payoff function be p : R � R D � R D !R þ , such that π(z • , z 0,0 , n 0,0 ) is the payoff to the focal individual with phenotype z • when the collection of average phenotypes among all other individuals is z 0,0 and the collection of environmental state variables across all patches is n 0,0 .We assume that the fitness of this focal individual can, in turn, be written as a function w : is a (D + 1)-dimensional vector collecting the payoff π • to the focal individual, the average payoff π 0 to a patch neighbour, and the average payoff π k to an individual from each patch k 6 ¼ 0.
As an argument to π 0 in Eq (16), we used z n 0;0 to denote a vector that is equal to z 0,0 except for its first entry, which is given by z n 0;0 ¼ z � =ðN À 1Þ þ ðN À 2Þz 0;0 =ðN À 1Þ, i.e. by the average trait among the neighbours of a neighbour of the focal individual.This captures the notion that the focal individual can influence the payoff of its neighbours.
Eq (15) allows individual fitness to depend on the payoff of all the individuals of its generation in an arbitrary way.This said, in most applications the survival and fecundity of an

Box 2. Multi-dimensional environment
We generalize s w (z) and s e (z) to the case where there are n e > 1 environmental state variables.We denote by n k;t ¼ ðn is the effect of environmental variable j in the focal patch on environmental variable i in patch k, t.The vector C(h) has i-th element given by C i (h), which is the Fourier transform of Under the infinite island model of dispersal, where R k,t = 0, e i,k,t = 0, and c i,0 j,k = 0 for all k 2 G except k = 0, Eq (II.D) reduces to eqs.15-16 of [41].
individual depend only on its own payoff.In this case, fitness may be written as where s : R þ !R þ and f : R þ !R þ are survival and fecundity functions, respectively, and quantities with R as a superscript are defined so that s Suppose we set survival s to zero in Eq (17).In that case, we obtain the fitness function of the classical Wright-Fisher process (e.g.eq. 3 in [50], in the absence of environmental effects and for a circular stepping-stone model).More generally, if s is a positive constant and payoffs only influence fecundity f, Eq (17) implements a form of "death-birth" updating protocol (i.e.where individuals sampled at random to die are replaced by the offspring of selected individuals according to payoff, e.g.[36,51]).Conversely, a "birth-death" updating is obtained by setting f to a positive constant and letting payoff influence survival s only.Eq (17) will constitute a useful platform to explore more specific examples later, even though many of our results hold for the more general relationship between payoff and fitness given by Eq (15).

Selection under payoff-mediated fitness.
If fitness is of the form of Eq (15), the selection gradient can be written as k k;0 @pðz � ; z 0;0 ; n 0;0 Þ @z k;0 |ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl } where and (see Appendix D in S1 Text for a derivation).To understand Eq (18), it is first useful to interpret λ k (Eq 20) as a coefficient of fitness interdependence through payoffs.Specifically, λ k measures the effect on the fitness of the focal individual of a change in the payoff of an individual at a distance k, relative to the effect of the payoff of the focal individual on its own fitness.When positive, λ k can thus be interpreted as the strength of competition, as it indicates how much an increase in the payoffs of an individual at distance k reduces the fitness of the focal individual.With this in mind, the coefficient κ k,t (Eq 19) can be thought of as a measure of relatedness scaled to competition (or scaled relatedness for short, e.g.[7,52,53]; with Eq 19 extending to isolation by distance the formalization of this concept found in [54]).To see this, suppose that the evolving trait helps current patch neighbours (i.e.such that @π(z • , z 0,0 , n 0,0 )/ @z k,0 > 0) and consider the numerator of κ 0,0 on the first line of Eq (19).This numerator consists in the relatedness R 0,0 of a focal individual towards current patch neighbours, discounted by an inclusive fitness effect through increased competition, which is due to the increase in neighbours' payoff caused by the focal individual helping them.The inclusive fitness effect consists of three relatedness weighted fitness effects: (i) −λ 0 /(N − 1), which captures the increase in competition experienced by the focal individual itself (where 1/(N − 1) is the frequency of the focal individual among the neighbours of an individual it helps); (ii) −λ 0 (N − 2) R 0,0 /(N − 1), which captures the increase in competition experienced by the relatives in the focal patch that are distinct from the individuals the focal individual is helping directly; and finally (iii) the total effect À P j2Gn0 λ j R j;0 , which captures the increase in competition experienced by relatives in other patches.The other two lines of Eq (19) are interpreted similarly: the relatedness towards neighbours from patch k, t discounted by an inclusive fitness effect on all recipients experiencing a change in fitness stemming from the focal changing the payoff of individuals in patch k, t (the denominator in Eq (19) standardizes these effects relative to the inclusive competitive effects induced by the focal individual changing its own payoff).More broadly, κ k,t can be interpreted as the number of units of its own payoff that the focal individual is willing to exchange with a unit of payoff accruing to a randomly sampled individual from patch k, t without changing the mutant's probability of fixation at a singular value z*.The scaled-relatedness coefficient κ k,t can thus be seen as the genetic value of other individuals in patch k, t relative to that of the focal individual in units of payoff.For (k, t) 6 ¼ (0, 0), κ k,t then measures the net social discount rate of temporally delayed and spatially extended fitness effects [43].
From the considerations above, Eq (18) can be read as an inclusive fitness effect at the payoff level.That is, selection depends on how the focal individual influences its own payoff and the payoff of all other individuals across patches, now and in the future, weighted by their scaled relatedness κ k,t .For recipients in the future (t � 1), payoff effects are mediated by how the focal individual perturbs the environment in each patch (via e k−j,t in Eq 18), and in turn by how such environmental perturbation influences payoffs (via @π/@n j,0 in Eq 18).
In Appendix E in S1 Text, we use Fourier analysis to compute scaled relatedness κ k,t for the fitness model given by Eq (17) and an arbitrary dispersal distribution m k .Our results are summarized in Box 3.For example, we obtain that under a Wright-Fisher process, holds, where is the probability that, under neutrality, a gene descending from the focal individual is in patch k at t > 0 steps in the future (which depends on the characteristic function MðhÞ of the dispersal distribution, Table 1).The collection of these probabilities, p t ¼ ðp k;t Þ k2G , can thus be seen as the distribution of a standard random walk on the set of patches with step distribution characterized by the dispersal probabilities m k .When such a random walk leads to a

Box 3. Scaled-relatedness coefficients
With fitness given by Eq (17), we show in Appendix E in S1 Text that where as usual, all functions are evaluated at the resident trait value z, a prime denotes a derivative, and L k ðXÞ is the inverse Fourier transform of a function X(h) at k 2 G (Eq I. D of Box 1).The functions F and G t are defined as For fecundity effects under a Wright-Fisher process (s ¼ s 0 ¼ 0), the above reduces to F(h) = 0 and G t ðhÞ ¼ f 0 MðhÞ t , which yields Eq (21) of the main text when f 0 ¼ 1 (i.e. the payoff is directly fecundity).
Using Eq (III.A), we also show in Appendix F in S1 Text that the summary statistic K, for selection on environmentally mediated social interactions under local interactions (Eq 24), is given by ; ðIII:CÞ where For fecundity effects under a Wright-Fisher process (s ¼ s 0 ¼ 0) the summary statistic K simplifies to MðhÞCðÀ hÞ 1 À CðÀ hÞMðhÞ ; ðIII:EÞ while for survival effects (f 0 ¼ 0), it simplifies to probability p k,t that is greater than under a uniform distribution (i.e. when p k,t > 1/D holds), Eq (21) indicates that scaled relatedness κ k,t is positive, so that selection favours environmental transformations that increase payoffs in patch k, t.Conversely, selection favours environmental transformations that decrease payoffs in patches where p k,t is smaller than under a uniform distribution (i.e. when p k,t < 1/D holds).Which patches are those depends on the dispersal distribution (compare Fig 2E with  The selection gradient (18) recovers previous results of social evolution theory, which provides a robustness check of our calculations.To see these connections, assume that there are no ecologically meditated interactions, that is, that @π(z • , z 0,0 , n 0,0 )/@n 0,0 = 0 holds, that −C(z) = @π(z • , z 0,0 , n 0,0 )/@z • < 0 is a net payoff cost to self, and that B k (z) = @π(z • , z 0,0 , n 0,0 )/@z k,0 is a payoff benefit to individuals at distance k, which is typical of models under the heading of the evolution of "cooperation" or "altruism".Further, suppose that individuals interact only with other individuals at a single distance k so that the selection gradient is proportional to −C(z) + κ k,0 B k (z).Then, Eq (18) entails that selection favours such a helping behaviour if κ k,0 > C(z)/ B k (z) holds, that is, if the scaled-relatedness coefficient is greater than the cost-to-benefit ratio.For a Wright-Fisher process, κ k,0 reduces to −1/(DN − 1) � 0 (Eq 21), which is always nonpositive.Thus, helping cannot spread regardless of population structure because κ k,0 > C(z)/ B k (z) cannot be satisfied as long as B k (z) > 0. This result was first derived for a lattice-structured population for D ! 1 and k = 0 in [55], and for finite D and any k under a circular one-dimensional habitat in [5] (chapter 8, and generalized to other abelian group structures in [44] and [43]).More generally, the scaled-relatedness coefficient given in Box 3 recovers established conditions for the spread of helping and harming behaviour in lattice-structured populations under different biological assumptions, such as for survival effects or for fecundity effects with overlapping generations (e.g.[15,36,51,56]; see [43] for the explicit connections to this literature).Finally, in the presence of ecologically mediated interactions, so that @π(z • , z 0,0 , n 0,0 )/@n 0,0 6 ¼ 0, Eq (18) recovers eq.(A.21) in [2] which holds for a Wright-Fisher process (to see this correspondence, set s k,t = N(@π(z • , z 0,0 , n 0,0 )/@n 0,0 )e k,t in Eq 18).

Local interactions.
In Eq (18), payoffs can depend on the traits expressed in and the environmental variables of all patches.In many instances, payoff effects can reasonably be assumed to be local, i.e. the payoff of an individual depends only on the traits and the environment of its patch.In this case, the payoff to the focal individual can be written as π(z • , z 0,0 , n 0,0 ) = π(z • , z 0,0 , n 0,0 ), and the selection gradient Eq (18) reduces to sðzÞ / @pðz � ; z 0;0 ; n 0;0 Þ @z � þ @pðz � ; z 0;0 ; n 0;0 Þ @z 0;0 k 0;0 |ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl } / s w ðzÞ þ @pðz � ; z 0;0 ; n 0;0 Þ @n 0;0 NK |ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl } / s e ðzÞ ; ð23Þ where summarizes selection on environmentally mediated social interactions (see Appendix F in S1 Text for a derivation and Eq III.C in Box 3 for more details).When K = 0, selection is blind to the effects of the trait on the environment, even if the environment affects payoffs (i.e. even if @π(z • , z 0,0 , n 0,0 )/@n 0,0 6 ¼ 0).When K > 0, selection favours trait values that improve the environment (i.e. the payoff in the future increases).Conversely, when K < 0, selection favours trait values that deteriorate the environment (i.e. the payoff in the future decreases).Which of these outcomes unfolds depends on the interaction between extended phenotypic effects e k,t and scaled-relatedness coefficients κ k,t , with K > 0 when e k,t and κ k,t tend to be of the same sign (i.e. when P 1 t¼1 P k2G e k;t k k;t > 0 holds), and K > 0 when they tend to be of opposite sign (i.e. when P 1 t¼1 P k2G e k;t k k;t < 0 holds).In the next section, we explore through an example how this interaction depends on dispersal and how environmental state variables of different patches influence one another.

Inter-temporal helping and harming through a lasting commons
To gain more specific insights into how isolation by distance influences the way selection shapes environmentally mediated interactions, consider a scenario where the environmental variable is some lasting commons (e.g. a common-pool resource or a toxic compound) that can move in space, and whose production depends on an evolving trait that is individually costly to express.We assume that the commons is a "good" when the environmental variable takes positive values (n > 0) and a "bad" when it takes negative values (n < 0).We also assume that the trait leads to the former when z > 0 and to the latter when z < 0, where the trait space is assumed to be Z ¼ ðÀ 1; 1Þ.The trait can thus be broadly thought of as environmentally mediated helping (increasing survival and reproduction to recipients) when z > 0, and as environmentally mediated harming (decreasing survival and reproduction to recipients) when z < 0.
We assume that fitness takes the form of Eq (17) with payoff given by where B > 0 and C > 0 are parameters that respectively tune the effects of the environmental variable in the focal patch n 0,0 and of the modifying trait z • of the focal individual on the payoff of the focal individual.These effects also depend on the shape parameters α B and α C , which we assume are positive integers, with α B odd (e.g.α B = 1) and α C even (e.g.α C = 2).Thus, the local commons increases (resp.decreases) payoffs when n 0,0 > 0 (resp.n 0,0 < 0) holds, but any trait expression, that is, any z • away from 0, is individually costly and reduces individual payoff, and thus fitness.We also assume that costs increase more steeply than benefits, i.e. that α C > α B holds.
Meanwhile, how the trait modifies the commons is determined by the environmental map g (Eq 2), which we assume is given by Eq (26) states that the commons changes from one demographic time point to the next due to three processes.First, the commons is modified or "produced" in a patch according to a function P : R !R of the average trait expressed in that patch, i.e. z R 0;0 ¼ z � =N þ ðN À 1Þz 0;0 =N in the focal patch and z k,0 otherwise (with k 6 ¼ 0).We assume that P(0) = 0 holds, and that P is monotonically increasing, i.e.P 0 (z) > 0 for all z 2 R, where P 0 (z) denotes the derivative of P with respect to z.Second, each unit of commons "diffuses" or moves with probability d k to a patch at a distance k from its source patch, where the probability distribution defined by d k can be thought of as the environmental equivalent of the dispersal probability distribution m k .We let DðhÞ denote the Fourier transform of this distribution for future use (i.e.DðhÞ is to d k what MðhÞ is to m k , Table 1).Third, a unit of commons decays with rate 0 < � � 1 from one time step to the next.Substituting Eq (26) into Eq (3) indicates that in a monomorphic population for z, the dynamics of the commons stabilises to which is positive when z > 0, negative when z < 0, and whose absolute value increases as the rate of decay � decreases, as expected.This equilibrium is unique because g is linear in P, and stable because � > 0 holds.With fitness of the form of Eq (17) and payoff Eq (25) only depending on local interactions (as in section 3.3.3),we can use Eqs ( 23) and ( 24) to characterise selection.With κ k,t as given in Box 3, all that remains to be computed are the extended phenotypic effects, e k,t .To do so, we start by substituting Eq (26) into Eqs ( 12) and ( 13), to obtain ψ k = P 0 (z)d k /N, and c k = (1 − �)d k , which in turn yield CðhÞ ¼ P 0 ðzÞDðhÞ=N and CðhÞ ¼ ð1 À �ÞDðhÞ.Substituting these expressions into Eq (14) then gives E t ðhÞ ¼ ð1 À �Þ tÀ 1 P 0 ðzÞDðhÞ t =N, which leads to for the extended phenotypic effect e k,t on patch k, t, where Equation Eq (28) can be understood as follows.By marginally changing its trait value, a focal individual produces P 0 (z)/N additional units of commons.Each such unit decays with time according to (1 − �) t−1 , and ends up in patch k, t according to q k,t (Eq 29), which can be interpreted as the probability that a non-decaying unit of the commons modified in the focal patch is located in patch k t generations in the future.In fact, the collection q t ¼ ðq k;t Þ k2G yields the distribution of a standard random walk on the set of patches with step distribution characterised by d k .Extended phenotypic effects thus depend critically on the way the commons moves in space as captured by d k (see Fig 4 for examples of e k,t in a 1D model).
In turn, how selection depends on extended phenotypic effects is found by substituting Eqs ( 25) and (28) into Eq (23).From this, we obtain can be thought of as the total expected genetic value in units of payoff of all the individuals in the future that are affected by a unit of the commons in the focal patch.Since the term multiplying O in Eq (30) is positive (recall that α B is odd and hence α B − 1 is even), the selection gradient increases with O, and the greater O is, the greater the z favoured by selection.In fact, the selection gradient reduces to s(0) / O at z = 0 (under our assumptions about parameters and the commons production function P).This shows that in a population where the trait is initially absent (so that individuals have no effect on the commons), selection favours environmental modifications leading to a common good (z > 0) when O > 0, or to a common bad (z < 0) when O < 0. Put differently, selection favours environmentally mediated inter-temporal helping when, in the eyes of the focal individual, the recipient of such help on average has positive genetic value in units of payoff, and conversely, inter-temporal harming when it has negative genetic value.
Further insights can be generated if we assume the simple functional form P(z) = P 0 z.In this case, we obtain The singular value z* (32) is convergence stable under our assumption that the cost of trait expression increases faster than the benefits (i.e.α C > α B ). Therefore, z* is the mode of the stationary phenotypic distribution of the trait substitution sequence (we compute this distribution for the Wright-Fisher model in eq.A-136 in S1 Text).From Eq (32) it is clear that the absolute value of z* increases with the benefit-to-cost ratio B/C, with α B /α C , and with the environmental effect of the trait P 0 .However, whether z* is positive or negative (and thus whether helping or harming evolves), ultimately depends on the sign of O, i.e. whether the expected genetic value O of a modification to the commons is positive or negative in payoff units.
The impact of species dispersal and commons movement on O can be understood most easily by assuming that payoff influences fecundity under a Wright-Fisher process (i.e.f 0 > 0 and s 0 ¼ s ¼ 0 in Box 3).In this case, O can be expressed as (see Appendix G.1 in S1 Text) where cov(p t , q t ) is the covariance between the distributions of the random walks of genes p t ¼ ðp k;t Þ k2G (Eq 22) and commons q t ¼ ðq k;t Þ k2G (Eq 29).This covariance is positive when there is a positive association between gene lineages and the commons these lineages modify.
In other words, O is positive and helping is favoured when an environmental modification owing to the expression of a gene is most likely to be experienced by individuals living in the future and carrying identical-by-descent copies of that gene.Conversely, O tends to be negative and harming is favoured when this environmental modification is less likely to be experienced by future carriers.While Eq (33) offers intuition on the biological conditions leading to positive or negative values of O, this quantity is more readily computed by noting that O = �KN/P 0 (z) holds, where K is defined in Eq (24), and by substituting eqs.(A-127) and (A-128) in S1 Text into Eq (III.E) in Box 3. Doing so, we obtain ) indicate that spite tends to be favoured by: (i) high levels of dispersal in the evolving species; (ii) high levels of movement of the commons; (iii) high environmental decay �; and (iv) significant differences in the dispersal distance of the species and of the commons (e.g. when individuals disperse short distances but the commons move far away from their original patch).These conditions lead to a negative association between gene lineages and the commons these lineages modify.Conversely, altruism tends to be favoured when dispersal and movement are weak, environmental decay is low, and the distributions of the species' dispersal and of the commons' movement are similar (Fig 5A and 5B, white region).In fact, under weak dispersal and movement (so that m 0 = 1 − m and d 0 = 1 − d with m and d close to zero), and regardless of the dispersal and movement distributions, we have which is always positive when m and d are sufficiently small (see Appendix G.2 in S1 Text for a derivation of Eq 35).Under these assumptions, an individual's lineage and the commons originating from its patch will be strongly and locally associated.
To see the effects of isolation by distance on social evolution, we compare in Fig 6A the singular strategy z* for a 2D lattice where the evolving species and the commons follow a binomial model for distance travelled (detailed in Appendix G.5 in S1 Text), with the singular strategy z* under the island model (where dispersal and commons movement are uniform).This comparison illustrates how isolation by distance allows for a wider range of evolved social ).These differences are down to the fact that isolation by distance allows for a wider range in the covariance among genetic and commons movement (recall Eq 33).We also compared the singular strategy z* found by substituting Eq (34) into Eq (32) with results from individual-based simulations.In these simulations, each offspring mutates with probability 10 −4 , in which case a normally distributed deviation with mean 0 and standard deviation 10 −2 is added to the parental trait value.The only difference between these simulations and our analytical scenario is thus that multiple alleles can segregate due to mutation (rather than just two alleles under a trait substitution sequence, see [7] for finite populations).Nevertheless, we find an excellent fit between the convergence stable z* and the mean population trait value in simulations (Fig 6B ).
The case where payoff influences survival rather than fecundity (s 0 > 0 and f 0 ¼ 0 in Eq 17) is illustrated in Fig 7 (the expression of O for this case in terms of characteristic dispersal functions can be found in Appendix G.4 in S1 Text).This analysis reveals that harming tends to be favoured when baseline survival s is low, especially when environmental decay is also low (Fig 7C).This is because, otherwise, an individual may harm itself in the future.But apart from this, selection is not fundamentally different when payoff influences survival rather than fecundity in this model (i.e.under a birth-death vs. death-birth process).

Discussion
Our analyses characterise in two main ways the selection gradient on a trait that impacts the deterministic dynamics of environmental state variables that can be abiotic or biotic, which in turn feed back on individual survival and reproduction under isolation by distance.
First, we showed how selection on a trait due to its environmental effects can be understood in terms of how a focal actor influences the fitness of all future individuals via a modification to the environmental state variables that these individuals experience (Eqs 11 and II.D for the case of multiple environmental variables).The relevant trait-driven environmental modifications are formalized as extended phenotypic effects that quantify how a trait change in an actor individual in the present affects the environmental state variables in all patches at all future times (the e k,t effects, Eq 14).While extended phenotypic effects are typically thought to benefit the actor or related contemporaries directly [57,58], these effects in our model are all indirect, carrying over in space and time, thus influencing the fitness of future carriers of the actor's trait when dispersal is limited.The associations between environmental and genetic variation that are necessary for selection to target trait-driven environmental modifications are given by the product between the extended environmental effects e k,t and the relatedness coefficients R k,t (see Eq 11), both of which can be efficiently computed using Fourier transforms (Eqs 14 and II.E and II.F for a multivariate environment).These gene-environment associations indicate that selection favours traits or behaviours with environmental effects such that, when expressed by a focal individual, the environmental effects increase (resp.decrease) the fitness Panel C: Parameter region where selection favours the evolution of helping (O > 0) or harming (O < 0) under survival effects with adult survival probability s on the x-axis and environmental decay � on the y-axis (O computed from Eq 31 using eq.A-129 in S1 Text; other parameters: same as in Fig 4B, i.e. under long-range movement of the commons).See S1 Data for how to generate these figures using Mathematica.https://doi.org/10.1371/journal.pcbi.1012071.g007 of future individuals that are more (resp.less) related to the focal than other individuals at that same future generation.
The second version of the selection gradient that we derived is based on the extra assumption that interactions between individuals are mediated through a payoff function (see Eq 17), as in most of traditional evolutionary game theory (e.g.[8,59,60]).Selection on a trait due to its environmental effects can still be viewed from an actor-centered perspective, but this time at the payoff rather than at the fitness level.Specifically, selection can be quantified in terms of how a focal individual influences the payoff of all future individuals via modifications to the environment these individuals experience, weighted by the relatedness between these individuals and the focal, now scaled to take competition into account (the term proportional to s e (z) in Eq 18, with scaled relatedness given in Eqs 19 and 20).The concept of scaled relatedness is useful because it summarizes in a single quantity (here one for each spatial and temporal distance) all the consequences of interactions among related individuals for indirect selection [7,52,53,61].That is, scaled relatedness balances, on the one hand, the positive effects of boosting the reproductive success of relatives in a particular spatial position, with, on the other hand, the negative effects of increasing competition for these relatives by affecting the reproduction and survival of others across the habitat.The increase of kin competition can be strong enough to offset the indirect benefits of social behaviour when social interactions occur among contemporaries (and generations do not overlap, e.g.[5,55]).Because the strength of kin competition depends on the specifics of the life cycle (such as whether generations overlap or not, or whether payoff influences fecundity or survival), the evolution of direct social interactions is sensitive to such assumptions (see [52] for a review).This is notably the case under isolation by distance, where the evolution of altruism crucially depends on whether reproduction is modelled as a "death-birth" or a "birth-death" process (e.g.[62][63][64]).The "death-birth" process is akin to iteroparous reproduction with fecundity effects, and the resulting life cycle can sustain altruism; the "birth-death" process is akin to iteroparous reproduction with survival effects, and the life cycle inhibits altruism via increased kin competition (e.g.[15,17,18,56]).
In contrast, we have found that in our model of environmentally mediated social interactions through a lasting commons, whether selection favours the evolution of helping or harming depends weakly on whether payoff influences survival or fecundity.There are two explanations for this.The first is that, because of environmental legacy, the effects on recipients are felt in the distant future, which decreases the competition among the focal's own offspring [41,65].The second explanation is that, in our model, individuals and their environmental effects can move in space independently, further dissociating the positive and negative effects of interactions among relatives.This decoupling between benefits and costs means that natural selection can readily favor either posthumous altruism or posthumous spite in our model with non-overlapping generations.Which of these behaviours evolves depends on whether the combination of dispersal pattern and commons movement cause environmental effects to fall predominantly on individuals that are more or less related than average in the future.
Our findings on environmentally mediated posthumous spite merit further discussion as existing models suggest that the conditions for the evolution of spite are more restrictive than those for altruism.By spite, we refer here to a trait or behaviour whose expression decreases the individual fitness of both its actor and recipients.This is a strong form of spite, reducing the actor's fitness (chapter 7 in [5]), which contrasts with the more commonly explored scenarios of weak spite (where the behaviour directly increases the actor's fitness, e.g.[66,67]).The evolution of strong spite typically relies on the existence of mechanisms by which individuals can evaluate their relatedness with social partners and thus behave according to some kin or type-recognition mechanism (e.g.[68][69][70][71][72]).By contrast, in our model spite is indiscriminate: an individual deteriorates the environment in the future without paying attention to recipients' identities, even if this comes at a cost in the present.In the absence of generational overlap (e.g.Wright-Fisher model), such spite is entirely posthumous.With this in mind, it is noteworthy that spite can evolve even when local populations are not small (e.g. of local size 50, Fig 5).More broadly, our results illustrate how environmentally mediated social interactions under isolation by distance can evolve to be as relevant for fitness as direct social interactions.
The two main assumptions of our model are that fitness and environmental effects are homogeneous in space and time, and that environmental dynamics are deterministic.These assumptions are common to previous models interested in environmental or ecological changes in space that are evolutionarily driven, and particularly to those where individuals produce an environmental commons that moves according to a diffusion process (e.g.[73][74][75][76][77]).These models further assume a separation of time scales between demography and the commons such that in between the reproduction, death, or dispersal of any individual in the entire population, the dynamics of the commons reaches a stable distribution across the landscape.By contrast, here we have assumed that the mutation process, rather than the demographic process, is slow compared to environmental dynamics.Reproduction, death, or dispersal can occur on a similar time scale than environmental dynamics in our model, as is usually the case in ecological models (e.g.[78,79]).As a result, even though environmental dynamics are described by a deterministic system (Eq 9, also as in most ecological models, [78,79]), realised environmental dynamics fluctuate randomly on a similar time scale than unavoidable genetic fluctuations owing to finite patch size.The next challenge would be to consider a fully stochastic system for the environmental variables (i.e. to extend Eq 9 to the dynamics of a probability distribution).This would be especially useful to investigate the effects of demographic stochasticity in response to trait evolution [25], and to model, for instance, environmentally mediated evolutionary suicide or rescue.Our framework may nevertheless provide a suitable approximation to cases of demographic and environmental stochasticity (with Eq 9 giving the expectation in state variable at the next time step, conditional on the states of at the previous step).This approach has been shown to work well under the island model of dispersal provided that patches were not too small and dispersal not too limited [41].It would be interesting to investigate how this holds up under isolation by distance.Finally, we considered the evolution of a single trait but our analyses apply directly to the case where multiple traits with environmental effects evolve jointly.To investigate such a case, one would first use our equations to obtain the selection gradient on each trait separately (applying Eq 5 to each trait under focus), and then jointly consider these gradients to perform a standard multidimensional analysis of directional selection (e.g.[80,81]).

Fig 1 .
Fig 1.A model for environmentally mediated social interactions in space and time.Schematic description of our model for a onedimensional lattice habitat (see sections 2.1-2.2 for details).Each patch k 2 {. .., D − 1, 0, 1, . ..} at time t 2 {0, 1, . ..} in the past is characterized by an environmental state variable n k,t (represented here by a cloud, e.g.water level, concentration of a pollutant, density of a resource), and the average trait value z k,t expressed by the individuals it carries (e.g.water absorption rate, detoxifying capacity, handling time; individuals represented here as palms).The environmental state n 0,0 of the focal patch k = 0 at time t = 0 depends on all environmental states and traits of the previous generation according to the environmental map g (blue dashed arrows, Eq 2).In turn, the fitness of a focal individual with trait z • (in yellow) depends on all environmental states and traits expressed in its own generation according to the fitness function w (orange arrows, Eq 1).The two functions g and w thus characterise how evolutionary and environmental dynamics interact with one another through dual inheritance of traits and environmental state variables.https://doi.org/10.1371/journal.pcbi.1012071.g001

Fig 2 .
Fig 2. Dispersal distribution, relatedness, and scaled relatedness in a 1D lattice model under short and long-range dispersal.Panels A-B: Dispersal distribution m k in a lattice-structured population in a one-dimensional habitat (with D 1 = 51).An offspring leaves its natal patch with probability 1 − m 0 = m = 0.8 and disperses to a patch at a Manhattan distance that follows a truncated binomial distribution (eq.A-7 in Appendix B.1 in S1 Text) with mean � λ m ¼ 1:9 in panel A, leading to short-range dispersal, and � λ m ¼ 15 in panel B, leading to long-range dispersal.The distance dispersed along each habitat dimension is uniformly distributed across all dimensions and directions (Appendix B.1 in S1 Text for details).Panels C-D: Relatedness R k,t for the dispersal distributions shown in panels A and B, respectively (using Eq 8 with patch size N = 20 and no adult survival s ¼ 0).Panel C highlights how relatedness decays in time and space, becoming negative away from the focal deme when dispersal is short-range.In contrast, in panel D, where dispersal is long-range, relatedness is negative at intermediate and large distances, thus leading to a multimodal distribution of relatedness values.Panels E-F: Scaled relatedness κ k,t for the dispersal distributions shown in panels A and B, respectively, under a Wright-Fisher model with fecundity effects (using Eq 21 with patch size N = 20).The trend of scaled relatedness is similar as that for relatedness.See S1 Data for how to generate these figures using Mathematica.https://doi.org/10.1371/journal.pcbi.1012071.g002

Fig 3 .
Fig 3. Dispersal distribution, relatedness, and scaled relatedness in a 2D lattice model under short and long-range dispersal.Panels A-B: Dispersal distribution m k in a two-dimensional habitat (with D 1 = D 2 = 13).An offspring leaves its natal patch with probability 1 − m 0 = m = 0.8 and disperses to a patch at a Manhattan distance that follows a truncated binomial distribution with mean � λ m ¼ 1:5 in panel A, leading to short-range dispersal, and � λ m ¼ 11 in panel B, leading to long-range dispersal (see Appendix B.2 in S1 Text for details).Panels C-D: Relatedness R k,0 from the dispersal distributions shown in panels A and B, respectively (using Eq 8 with patch size N = 20 and no adult survival s ¼ 0).Panels E-F: Scaled relatedness, κ k,10 in panel E and κ k,1 in panel F, from the dispersal distributions shown in panels A and B, respectively, for a Wright-Fisher model (using Eq 21 with patch size N = 20).See S1 Data for how to generate these figures using Mathematica.https://doi.org/10.1371/journal.pcbi.1012071.g003 the payoff to self, the average payoff to a patch neighbour, and the average payoff to an individual from each patch other than the focal, i.e. as wðz � ; z 0;0 ; n 0;0 Þ ¼ wðπðz � ; z 0;0 ; n 0;0 ÞÞ; ð15Þ where πðz � ; z 0;0 ; n 0;0 Þ ¼ ðpðz � ; z 0;0 ; n 0;0 Þ |ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl } p � ; pðz 0;0 ; z n 0;0 ; n 0;0 Þ |ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl } p 0 ; . . .; pðz k;0 ; z R k;0 ; n k;0 Þ |ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl } p ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl } reproduction into spots left open by deaths ; ð17Þ

e
kÀ j;t @pðz � ; z 0;0 ; n 0;0 Þ @n j;0 Nk k;t |ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl } / s e ðzÞ ; ð18Þ Fig 2F for short and long-range dispersal in a 1D lattice, and Fig 3E with Fig 3F for short and long-range dispersal in a 2D lattice).

Fig 4 .
Fig 4. Extended phenotypic effects in a 1D lattice model under short and long-range movement of the commons.Panel A: When the commons moves locally, extended phenotypic effects e k,t decay in time and space away from the focal deme (from Eq 28 with D 1 = 31, movement probability d = 0.6 and expected movement distance � λ d ¼ 1:55, see Appendix G.5 in S1 Text for details on how movement is modelled; production function P(z) = Nz, i.e. each unit of z contributes to one unit of resource; decay rate � = 0.2; other parameters: N = 20, m = 0.3, � λ m ¼ 1:55).Panel B: In contrast, when the resource moves at greater distances, extended phenotypic effects e k,t are greatest further away from the focal deme (from Eq 28 with movement parameters d = 0.98 and � λ d ¼ 8; production function P(z) = Nz; decay rate � = 0.5; other parameters: same as Fig 2A).See S1 Data for how to generate these figures using Mathematica.https://doi.org/10.1371/journal.pcbi.1012071.g004

Fig 5 .
Fig 5. Selection favours altruism or spite depending on dispersal of the evolving species and the commons in a 2D lattice model.Panels A-B: Regions of dispersal parameters leading to the evolution of altruism, O > 0 (in white), or of spite, O < 0 (in gray) for an example in a 2D lattice model (with D 1 = D 2 = 13 and N = 50) under a Wright-Fisher life-cycle with fecundity effects (with O computed from Eq 34).Panel A: Combination of dispersal probability of the evolving species m = 1 − m 0 (x-axis) and of the commons d = 1 − d 0 (y-axis) for different levels of environmental decay � in different shades of gray (� = 0.1, 0.5, 1) with expected dispersal distance fixed ( � λ m ¼ 1:54 and � λ d ¼ 8).This shows that spite is favoured by high levels of dispersal and environmental decay.Panel B: Combination of expected dispersal distance of the evolving species � λ m (x-axis) and of the commons � λ d (y-axis) for different levels of environmental decay � in different shades of gray (see panel A for the legend) with dispersal probability fixed (m = 0.98 and d = 1).This shows that spite is favoured by dispersal asymmetry between the evolving species and the commons.Panels C-D: Evolution of spite in individual-based simulations under a Wright-Fisher life-cycle with fecundity effects (with D 1 = D 2 = 13, N = 50, m = 0.3, � λ m ¼ 1:54, d = 1, � λ d ¼ 8, B = 2, α B = 1, C = 1, α C = 4, P(z) = Nz; for mutation, the trait mutates during reproduction with probability 10 −4 , in which case a normally distributed deviation with mean 0 and standard deviation 10 −2 is added to the parental trait value).Panel C shows the average trait z in the population; panel D shows the average commons level or environmental variable n (with simulations in full lines-see S1 Code-and analytical prediction in dashed lines-from Eq 32 for z and 27 for n).See S1 Data for how to generate these figures using Mathematica.https://doi.org/10.1371/journal.pcbi.1012071.g005

Fig
Fig 5A and 5Bgive the sign of O under a binomial model for the distance of both the dispersal of the focal species and the movement of the commons (see Appendix G.5 in S1 Text for details).These figures show that such model of dispersal allows for both positive and negative values of O and thus for the evolution of both inter-temporal helping and harming.Here, helping corresponds to posthumous altruism and harming to posthumous spite since an individual can never obtain direct benefits from its own trait effects on the environment, and these effects are only felt by the next generation (as generations are non-overlapping under a Wright-Fisher process).More generally, numerical explorations of O (Fig 5Aand 5B) indicate that spite tends to be favoured by: (i) high levels of dispersal in the evolving species; (ii) high levels of movement of the commons; (iii) high environmental decay �; and (iv) significant differences in the dispersal distance of the species and of the commons (e.g. when individuals disperse short distances but the commons move far away from their original patch).These conditions lead to a negative association between gene lineages and the commons these lineages modify.Conversely, altruism tends to be favoured when dispersal and movement are weak, environmental decay is low, and the distributions of the species' dispersal and of the commons' movement are similar (Fig 5Aand 5B, white region).In fact, under weak dispersal and movement (so that m 0 = 1 − m and d 0 = 1 − d with m and d close to zero), and regardless of the dispersal and movement distributions, we have

Fig 6 .Fig 7 .
Fig 6. Isolation by distance allows for the evolution of a wider range of social behaviours than the island model under environmental feedback.Panel A: Singular trait value z* from Eq 32 for m = 0.75 (in black) and m = 0.2 (in gray).Other parameters: D 1 = D 2 = 13, so that D = 169, N = 50, � λ m ¼ 1:54, d = 0.99, B = 2, α B = 1, C = 1, α C = 4, P(z) = Nz.Dashed lines show the singular trait value with the same parameters but under the island model of dispersal for both the species and the commons (from Eq 32 with O given by eq.(A-126 in S1 Text) derived in Appendix G.3). Panel B: Observed vs. predicted equilibrium trait value in individual-based simulations running for 20 000 generations under different expected dispersal distance of the commons � λ d leading to altruism (z > 0) and spite (z < 0).Other parameters: same as in Fig 5C and 5D.The prediction is shown as a dashed line (from Eq 32) with grey region around for twice the standard deviation obtained from the stationary phenotypic distribution (from eq.A-136 in S1 Text).Observed values of the trait average in the population are shown as black dots for the average from generation 5 000 to 20 000, with error bars for standard deviation over the same 15 000 generations.Simulations were initialised at the predicted convergence stable trait value.S1 Code for simulation code and S2 Data for simulation results.https://doi.org/10.1371/journal.pcbi.1012071.g006

Table 1 . Key general symbols.
� w k ðhÞ conjugate of character function (Eq I.E in Box 1) L k ðF Þ inverse transform of F at k (Eq I.D in Box 1) e k,t extended phenotypic effect of the focal individual on the environmental state variable in patch k at time t in the future (Eqs 10 and 14) ψ k focal individual's effect on the environmental state variable of patch k over one generation (Eq 12) C(h) Fourier transform (Eq I.B in Box 1) of ψ k c k effect of the environmental state variable of one patch on the environmental state variable of another patch at distance k over one generation (Eq 13) CðhÞ Fourier transform (Eq I.B in Box 1) of c k π individual payoff function (section 3.3.1)λk coefficient of fitness interdependence through payoffs (Eq 20) κ k,t scaled-relatedness coefficient between the focal individual and another individual from patch k, t (Eq 19)p k,t probability that a gene descending from the focal individual is in patch k at t > 0 steps in the future (Eq 22) Key symbols from example section 3.4 B, α B parameters tuning the effects of the environmental variable on the payoff of the focal individual (Eq 25) C, α C parameters tuning the effects of the trait z • of the focal individual on the payoff of the focal individual (Eq 25) P commons production function (Eq 26) � commons decay rate (Eq 26) d k probability that a unit of commons moves to a patch at a distance k from its source patch (Eq 26) DðhÞ Fourier transform (Eq I.B in Box 1) of d k q k,t probability that a non-decaying unit of the commons modified in the focal patch is located in patch k, t generations in the future (Eq 29) O expected genetic value in units of payoffs of all the individuals in the future that are affected by a unit of the commons in the focal patch (Eq 31) https://doi.org/10.1371/journal.pcbi.1012071.t001

for a summary dia- gram of our model). We denote by z • the realized value of the trait of the focal individual, and by z k,t the realized average trait of individuals other than the focal individual living in patch k, t
|ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl } 1;k;t ; n 2;k;t ; . . .n n e ;k;t Þ the vector of such variables in patch k, t (where n i;k;t 2 R is the value of the ith environmental variable).The fitness of the focal individual is now given by i ðz; n Þ for i ¼ 1; 2; ...; n e ; ðII:CÞ where n ¼ ðn; ...; nÞ is a vector of dimension D whose entries are all given by n ¼ ðn 1 ; ...; nn e Þ.With fitness given by Eq (II.A), selection on intra-temporal effects s w (z) remains unchanged, given by Eq (6a) with wðz � ; z 0;0 ; ñ0;0 ðz H ÞÞ substituted for w(z • , z 0,0 , n 0,0 (z H )). For selection on inter-temporal effects s e (z), carrying out mutadis mutandis the same calculations as we have for the one-dimensional case, yieldsse ðzÞ ¼ k2G e i;jÀ k;t @wðz � ; z 0;0 ; ñ0;0 Þ @n i;j;0 NR k;t ; ðII:DÞ where e i,k,t is the extended phenotypic effect on environmental variable i in patch k, t. h2G E i;t ðhÞ� w k ðhÞ; ðII:EÞ where E i;t ðhÞ is the i-th element of the vector Ẽ t ðhÞ ¼ ðE 1;t ðhÞ; E 2;t ðhÞ; . . .; E n e ;t ðhÞÞ, which is obtained from Ẽ t ðhÞ ¼ CðhÞ tÀ 1 CðhÞ: ðII:FÞ Here, the community matrix C(h) has its ij-th element given by C ij ðθÞ ¼ P k2G c i;0 j;k w k ðhÞ, where